using System;
using System.Drawing;
using System.Data;
using System.Data.Odbc;
using System.Collections;

namespace StaticUtility
{	
	/*
	public class PatchStats
	{
		PointF CircleCenter;
		double CircleRadius;
		int Nparticles;
		double[] Mean;
		double[,] Cov;
		double[] EigVals;
		double[,] EigVecs;
		int Dim;
		public PatchStats() {}
		public PatchStats(int Dim)
		{
			this.Dim = Dim;
			Mean = new double[Dim];
			Cov = new double[Dim,Dim];
			EigVals = new double[Dim];
			EigVecs = new double[Dim,Dim];
		}
	}
	*/

	/// <summary>
	/// Static utilities.
	/// </summary>
	public class Util
	{
		public Util(){}
		public static readonly double FourThirdsPi = 4.0*Math.PI/3.0;
		public static readonly double TwoThirdsPi = 2.0*Math.PI/3.0;
		public static readonly double HalfPi = Math.PI/2.0;
		public static readonly double ThirdPi = Math.PI/3.0;
		public static readonly double QuarterPi = Math.PI/4.0;
		public static readonly double TwoPi = Math.PI*2.0;
		public static readonly double Sqrt2 = Math.Sqrt(2.0);
		public static readonly float PiFloat = (float)Math.PI;
		public static readonly float TwoPiFloat = (float)Util.TwoPi;
		public static Random Rand = new Random();
		/// <summary>
		/// Creates a random vector of integers from smallestValue to smallestValue+length-1.
		/// </summary>
		/// <param name="length">The lenght of the random vector.</param>
		/// <param name="smallestValue">The smallest value in the random vector.</param>
		/// <param name="output">Returns a random vector of integers from smallestValue to smallestValue+length-1.</param>>
		public static void RandVec(int length, int smallestValue, int[] output)
		{
			int[] set0 = new int[length];
			for(int i=0; i<length; i++) set0[i] = smallestValue + i;

			int limit = length - 1;
			int selection;
			for(int i=0; i<length; i++)
			{
				selection = (int)Math.Floor(Rand.NextDouble()*(1.0+(double)(limit)));
				output[i] = set0[selection];
				set0[selection] = set0[limit--];
			}
		}
		/// <summary>
		/// Returns a fully saturated color from a location on the color-wheel from blue (0) to red(1).
		/// </summary>
		/// <param name="p">A value between 0 (min, blue) and 1 (max, red).</param>
		/// <returns></returns>
		public static Color BlueGreenYellowRed(double p)
		{
			return Color.FromArgb(
				Math.Max(0,Math.Min(255,(int)(255.0*(2.0-6.0*Math.Abs(p-0.80))))), 
				Math.Max(0,Math.Min(255,(int)(255.0*(2.0-6.0*Math.Abs(p-0.50))))), 
				Math.Max(0,Math.Min(255,(int)(255.0*(2.0-6.0*Math.Abs(p-0.20))))) );
		}
		private static double _c, _s;
		/// <summary>
		/// Returns a fully saturated color from a location on the color-wheel.
		/// </summary>
		/// <param name="p">An angle on the color wheel in radians.</param>
		/// <returns></returns>
		public static Color ColorWheel(double p)
		{
			_c = Math.Cos(p);
			_s = Math.Sin(p);
			if(_s>0.0)
			{
				if(_c>0.0)	//	pos sin, pos cos (red to yellow)
				{
					return Color.FromArgb( 
						(int)(255.0*Math.Min(1.0,(_c+_s))),
						(int)(255.0*_s), 
						0);
				}
				else		//	pos sin, neg cos (yellow to green)
				{
					return Color.FromArgb( 
						(int)(255.0*_s),
						(int)(255.0*Math.Min(1.0,(_s-_c))), 
						0);
				}
			}
			else
			{
				if(_c>0.0)	//	neg sin, pos cos (blue to purple)
				{
					return Color.FromArgb( 
						(int)(_c*255.0),
						0, 
						(int)(-_s*255.0));
				}
				else		//	neg sin, neg cos (green to blue)
				{
					return Color.FromArgb( 
						0,
						(int)(-_c*255.0), 
						(int)(-_s*255.0));
				}
			}
		}
		public static void ImprintHorizontalColorWheel(Bitmap bmp, RectangleF rectf)
		{
			Graphics g = Graphics.FromImage(bmp);
			float w = rectf.Width/100.0f;
			for(float i=0.0f; i<99.5f; i++)
				g.FillRectangle(new SolidBrush(ColorWheel(0.0628f*i)), rectf.Left+i*w, rectf.Top, w, rectf.Height);
			g.DrawRectangle(Pens.Black, rectf.X, rectf.Y, rectf.Width, rectf.Height);
			g.Dispose();
		}
		/// <summary>
		/// The color of a standard normally distributed 'z' value.
		/// </summary>
		/// <param name="z"></param>
		/// <returns></returns>
		public static Color BlueGreenGrayYellowRed(float z)
		{
			if(z>=1.0)	//	Arc from 1.0 (yellow) to red (>1.0)
			{
				return ColorWheel(Math.Max(0.0,Util.HalfPi-(z-1.0)));
			}
			else if(z<=-1.0)	//	Arc from -1.0 (green) to purple(<-1.0)
			{
				return ColorWheel(Math.Min(5.0,Math.PI-(z+1.0)));
			}
			else if(z>0.0)	//	Arc from white (0.0) to yellow (1.0)
			{
				_col = ColorWheel(Util.HalfPi);
				_white = (1.0-z)/2.0;
				return Color.FromArgb(
					_col.R + (int)((double)(255-_col.R)*_white),
					_col.G + (int)((double)(255-_col.G)*_white),
					_col.B + (int)((double)(255-_col.B)*_white) );
			}
			else	//	Arc from green (-1.0) to gray (0.0)
			{
				_col = ColorWheel(Math.PI);
				_white = (1.0+z)/2.0;
				return Color.FromArgb(
					_col.R + (int)((double)(255-_col.R)*_white),
					_col.G + (int)((double)(255-_col.G)*_white),
					_col.B + (int)((double)(255-_col.B)*_white) );
			}
		}
		public static double _white;
		public static Color _col;
		/// <summary>
		/// The color of a standard normally distributed 'z' value.
		/// </summary>
		/// <param name="z"></param>
		/// <returns></returns>
		public static Color BlueGreenWhiteYellowRed(double z)
		{
			if(z>=1.0)	//	Arc from 1.0 (yellow) to red (>1.0)
			{
				return ColorWheel(Math.Max(0.0,Util.HalfPi-(z-1.0)));
			}
			else if(z<=-1.0)	//	Arc from -1.0 (green) to purple(<-1.0)
			{
				return ColorWheel(Math.Min(5.0,Math.PI-(z+1.0)));
			}
			else if(z>0.0)	//	Arc from white (0.0) to yellow (1.0)
			{
				_col = ColorWheel(Util.HalfPi);
				_white = 1.0-z;
				return Color.FromArgb(
					_col.R + (int)((double)(255-_col.R)*_white),
					_col.G + (int)((double)(255-_col.G)*_white),
					_col.B + (int)((double)(255-_col.B)*_white) );
			}
			else	//	Arc from green (-1.0) to white (0.0)
			{
				_col = ColorWheel(Math.PI);
				_white = 1.0+z;
				return Color.FromArgb(
					_col.R + (int)((double)(255-_col.R)*_white),
					_col.G + (int)((double)(255-_col.G)*_white),
					_col.B + (int)((double)(255-_col.B)*_white) );
			}
			/*
			if(z>1.0f)
			{
				//	arc between red and yellow.
				z = (z-1.0f)*0.59f;	//	max at 1.5708 (pi/2)
				if(z>Util.HalfPi)
					z=Util.HalfPi;
					
				//	Weight red with sin(z) and yellow with cos(z)
				_s = Math.Sin(z);
				_c = Math.Cos(z);

				return Color.FromArgb( (int)(Math.Min(255.0,(_c+_s)*200.0)), (int)(_c*255.0), 0 );
			}
			//	Yellow to white
			else if(z>=0.0)
			{
				_white = (1.0-z)*250.0;
				//	linear transition between yellow to white.
				return Color.FromArgb( (int)(z*125.0+white), (int)(1.1*z*125.0f+white), (int)white );
			}
			else if(z>-1.0f)
			{
				z=-z;
				double white = (1.0-z)*250.0;
				//	linear transition between white to green.
				return Color.FromArgb( (int)white, (int)(z*125.0f+white), (int)white );
			}
			else //	 z < -1.0f
			{
				//	arc between blue and green.
				z = -(z+1.0f)*0.59f;	//	max at 1.5708 (pi/2)
				if(z>1.5708f)
					z=1.5708f;
					
				//	Weight blue with sin(z) and green with cos(z)
				double s = Math.Sin((double)z);
				double c = Math.Cos((double)z);

				return Color.FromArgb( 0, (int)(c*255.0), (int)(s*255.0) );
			}
			*/
		}
		public static void ImprintHorizontalBgwyr(Bitmap bmp, RectangleF rectf)
		{
			Graphics g = Graphics.FromImage(bmp);
			float w = rectf.Width/100.0f;
			for(float i=0.0f; i<99.5f; i++)
				g.FillRectangle(new SolidBrush(BlueGreenWhiteYellowRed((i-50.0f)/17.5f)), rectf.Left+i*w, rectf.Top, w, rectf.Height);
			g.DrawRectangle(Pens.Black, rectf.X, rectf.Y, rectf.Width, rectf.Height);
			g.Dispose();
		}
        public static void ImprintVerticalBgwyr(Bitmap bmp,RectangleF rectf)
        {
            Graphics g = Graphics.FromImage(bmp);
            float h = rectf.Height/100.0f;
            for(float i=0.0f;i<99.5f;i++)
                g.FillRectangle(new SolidBrush(BlueGreenWhiteYellowRed((i-50.0f)/17.5f)),rectf.Left,rectf.Top+i*h,rectf.Width,h);
            g.DrawRectangle(Pens.Black,rectf.X,rectf.Y,rectf.Width,rectf.Height);
            g.Dispose();
        }
		public static double Logistic(double z)
		{
			return 1.0/(1.0+Math.Exp(-z));
		}
		public static double InvLogistic(double p)
		{
			p = p*0.999998+0.000001;
			return -Math.Log( (1.0-p)/p );
		}
		/// <summary>
		/// Outputs a string that represents the unsigned integer "num" with "padding" digits (usually led by zeros).
		/// </summary>
		/// <param name="num"></param>
		/// <param name="padding"></param>
		/// <returns></returns>
		public static string PaddedIntegerString(uint num, uint padding)
		{
			if(num >= (uint)Math.Pow(10.0,(double)padding+1.0))
				return Convert.ToString(num);
			
			string output = "";
			
			for(uint i=0; i<padding; i++)
			{
				output = Convert.ToString(num%10) + output;
				num = num/10u;
			}
			return output;
		}
		/// <summary>
		/// For a given filename with the path, this removes all the information regarding the path (if the path info exists).
		/// </summary>
		/// <param name="filepath"></param>
		/// <returns></returns>
		public static string RemovePathInfo(string filepath)
		{
			int ct = filepath.Length-2;
			while(ct>=0 && filepath[ct]!='/' && filepath[ct]!='\\')
				ct--;
			return filepath.Substring(ct+1);
		}
		/// <summary>
		/// For a given filename with the path, removes the filename and the last '/' or '\\'
		/// </summary>
		/// <param name="filepath"></param>
		/// <returns></returns>
		public static string RemoveFilenameFromPath(string filepath)
		{
			int ct = filepath.Length-2;
			while(ct>=0 && filepath[ct]!='/' && filepath[ct]!='\\')
				ct--;
			if(ct>0)
				return filepath.Substring(0,ct);
			else
				return filepath;
		}
		/// <summary>
		/// For a given filename with the path, removes the extension.
		/// </summary>
		/// <param name="filepath"></param>
		/// <returns></returns>
		public static string RemoveExtension(string filepath)
		{
			int ct = filepath.Length-2;
			while(ct>=0 && filepath[ct]!='.')
				ct--;
			if(ct>0)
				return filepath.Substring(0,ct);
			else
				return filepath;
		}
		/// <summary>
		/// This configures the attributes of col to be compatible with the sqlTypeString.
		/// Unsupported SQL types:
		///		DECIMAL, DEC, NUMERIC, FIXED, DATE, TIME, TIMESTAMP, YEAR, SET
		/// </summary>
		/// <param name="col">The data column.</param>
		/// <param name="sqlTypeString">The string containing the SQL type definition.</param>
		public static void SetDataColumnType(DataColumn col, string sqlTypeString)
		{
			string colTypeString = sqlTypeString.ToUpper();
			if( colTypeString.IndexOf("DATETIME")>=0 )
			{
				col.DataType = typeof(DateTime);
			}
			else if( colTypeString.IndexOf("BIT")>=0 || colTypeString.IndexOf("BOOL")>=0 )	//	implicitly includes BOOLEAN
			{
				col.DataType = typeof(bool);
			}
			else if( colTypeString.IndexOf("BIGINT")>=0 )
			{
				if( colTypeString.IndexOf("UNSIGNED")>=0 )
					col.DataType = typeof(ulong);
				else
					col.DataType = typeof(long);
			}
			else if( colTypeString.IndexOf("TINYINT")>=0 )
			{
				col.DataType=typeof(byte);
			}
			else if(colTypeString.IndexOf("SMALLINT")>=0 )
			{
				if( colTypeString.IndexOf("UNSIGNED")>=0 )
					col.DataType = typeof(ushort);
				else
					col.DataType = typeof(short);
			}
			else if( colTypeString.IndexOf("INT")>=0)		//	implicitly includes MEDIUMINT (24-bit), and INTEGER
			{
				if( colTypeString.IndexOf("UNSIGNED")>=0 )
					col.DataType = typeof(uint);
				else
					col.DataType = typeof(int);
			}
			else if( colTypeString.IndexOf("ENUM")>=0 )
			{
				col.DataType = typeof(int);
			}
			else if( colTypeString.IndexOf("FLOAT")>=0 )
			{
				col.DataType = typeof(float);
			}
			else if( colTypeString.IndexOf("REAL")>=0 || colTypeString.IndexOf("DOUBLE")>=0 )
			{
				col.DataType = typeof(float);
			}
			else if( colTypeString.IndexOf("CHAR")>=0 || colTypeString.IndexOf("TEXT")>=0)	//	implicitly includes VARCHAR, TINYTEXT, MEDIUMTEXT, LONGTEXT 
			{
				col.DataType = typeof(string);
			}
			else if( colTypeString.IndexOf("BLOB")>=0 )	//	implicitly includes TINYBLOB, MEDIUMBLOB, LONGBLOB
			{
				col.DataType = typeof(byte[]);
			}
			col.AllowDBNull = colTypeString.IndexOf("NOT NULL")<0;
		}
		/// <summary>
		/// Returns a string of n comma separated question marks.
		/// </summary>
		/// <param name="n">The number of question marks.</param>
		/// <returns></returns>
		public static string CommaSepQuestionMarks(int n)
		{
			if(n<=0)
				return "";
			else
			{
				char[] output = new char[2*n-1];
				output[0] = '?';
				int i,ii;
				for(i=1; i<n; i++)
				{
					ii = 2*i-1;
					output[ii] = ',';
					output[ii+1] = '?';
				}
				return new string(output);
			}
		}
		/// <summary>
		/// Creates a clone fo the OdbcParameter.
		/// </summary>
		/// <param name="par">OdbcParameter  The clonee.</param>
		/// <returns>odbcParameter  The cloned.</returns>
		public static OdbcParameter OdbcParameterCopy(OdbcParameter par)
		{
			OdbcParameter output = new OdbcParameter(par.ParameterName, par.OdbcType, par.Size, par.SourceColumn);
			output.IsNullable = par.IsNullable;
			output.Direction = par.Direction;
			output.SourceVersion = par.SourceVersion;
			output.Precision = par.Precision;
			return output;
		}
		/// <summary>
		/// Converts a matrix to a string.
		/// </summary>
		/// <param name="mat"></param>
		/// <returns></returns>
		public static string MatrixToString(double[,] mat)
		{
			string output = "";
			int rows = mat.GetLength(0);
			int cols = mat.GetLength(1);	
			for(int i=0; i<rows; i++)
			{
				for(int j=0; j<cols; j++)
					output = output + "\t" + Convert.ToString(mat[i,j]);
				if(i<rows-1)
					output = output + "\n";
			}
			return output;
		}
		/// <summary>
		/// Converts a matrix to a string.
		/// </summary>
		/// <param name="mat"></param>
		/// <returns></returns>
		public static string MatrixToString(int[,] mat)
		{
			string output = "";
			int rows = mat.GetLength(0);
			int cols = mat.GetLength(1);	
			for(int i=0; i<rows; i++)
			{
				for(int j=0; j<cols; j++)
					output = output + "\t" + Convert.ToString(mat[i,j]);
				if(i<rows-1)
					output = output + "\n";
			}
			return output;
		}
		/// <summary>
		/// Computes the angle of the segment a-->b.
		/// </summary>
		/// <param name="a"></param>
		/// <param name="b"></param>
		/// <returns>A value in the range [0, 2pi]</returns>
		public static double Angle(PointF a, PointF b)
		{
			double output =  Math.Atan2( (double)(b.Y-a.Y), (double)(b.X-a.X) );
			if(output<0.0)
				output+=Util.TwoPi;
			return output;
		}
		/// <summary>
		/// angle1 - angle2, in radians
		/// </summary>
		/// <param name="angle1"></param>
		/// <param name="angle2"></param>
		/// <returns></returns>
		public static double AngleDistance(double angle1, double angle2)
		{
			double output = angle1-angle2;
			if(output > Math.PI)
				output -= TwoPi;
			else if(output < -Math.PI)
				output += TwoPi;
			return output;
		}
		private static int _i, _j, _iLen, _jLen;
		/// <summary>
		/// Matrix multiplication. output = mat * colVec
		/// </summary>
		/// <param name="mat">A matrix with as many columns as the number of rows in colVec.</param>
		/// <param name="colVec">A vector with as many rows as the columns of mat</param>
		/// <param name="output">A column vector the same number of rows as mat.</param>
		public static void MM(float[,] mat, float[] colVec, float[] output)
		{
			_iLen = colVec.Length;
			_jLen = mat.GetLength(0);
			for(_j=0; _j<_jLen; _j++)
			{
				output[_j]=0.0f;
				for(_i=0; _i<_iLen; _i++)
				{
					output[_j]+=mat[_j,_i]*colVec[_i];
				}
			}
		}
		/// <summary>
		/// Matrix multiplication. output = rowVec* mat
		/// </summary>
		/// <param name="rowVec">A vector with as many columns as the rows of mat</param>
		/// <param name="mat">A matrix with as many rows as the columns in rowVec.</param>
		/// <param name="output">A row vector the same number of columns as mat.</param>
		public static void MM(float[] rowVec, float[,] mat, float[] output)
		{
			_iLen = mat.GetLength(1);
			_jLen = rowVec.Length;
			for(_i=0; _i<_iLen; _i++)
			{
				output[_i]=0.0f;
				for(_j=0; _j<_jLen; _j++)
				{
					output[_i]+=rowVec[_j]*mat[_j,_i];
				}
			}
		}

		/// <summary>
		/// Matrix multiplication. output = mat * colVec
		/// </summary>
		/// <param name="mat">A matrix with as many columns as the number of rows in colVec.</param>
		/// <param name="colVec">A vector with as many rows as the columns of mat</param>
		/// <param name="output">A column vector the same number of rows as mat.</param>
		public static void MM(double[,] mat, double[] colVec, double[] output)
		{
			_iLen = colVec.Length;
			_jLen = mat.GetLength(0);
			for(_j=0; _j<_jLen; _j++)
			{
				output[_j]=0.0f;
				for(_i=0; _i<_iLen; _i++)
				{
					output[_j]+=mat[_j,_i]*colVec[_i];
				}
			}
		}
		/// <summary>
		/// Matrix multiplication. output = rowVec* mat
		/// </summary>
		/// <param name="rowVec">A vector with as many columns as the rows of mat</param>
		/// <param name="mat">A matrix with as many rows as the columns in rowVec.</param>
		/// <param name="output">A row vector the same number of columns as mat.</param>
		public static void MM(double[] rowVec, double[,] mat, double[] output)
		{
			_iLen = mat.GetLength(1);
			_jLen = rowVec.Length;
			for(_i=0; _i<_iLen; _i++)
			{
				output[_i]=0.0f;
				for(_j=0; _j<_jLen; _j++)
				{
					output[_i]+=rowVec[_j]*mat[_j,_i];
				}
			}
		}
		public class Gaussian
		{
			private static double absz, density, t, cumulative;
			public static double Density(double z)
			{
				return .39894228040143*Math.Exp(-0.5*Math.Pow(z,2.0));
			}
			public static double Probability(double z)
			{
				absz = Math.Abs(z);
				density = .39894228*Math.Exp(-.5*Math.Pow((absz),2.0));
				t = 1.0/(1.0+0.33267*absz);
				cumulative =
					1.0 - density*(0.4361836*t-0.1201676*Math.Pow(t,2.0)+0.937298*Math.Pow(t,3.0));

				if(z<0)
					return 1.0-cumulative;
				else
					return cumulative;
			}
		}
		/// <summary>
		/// Differences f into df using forward and backward values of f (except for first and last values of f).
		/// </summary>
		/// <param name="f"></param>
		/// <param name="df"></param>
		/// <param name="len"></param>
		public static void Diff(double[] f, double[] df, int len)
		{
			if(len>1)
			{
				df[0]=f[1]-f[0];
				df[len-1]=f[len-2]-f[len-1];
				for(int i=1; i<len-2; i++)
					df[i]=(f[i+1]-f[i-1])/2.0;
			}
		}
		/// <summary>
		/// Computes the orientation for each element of 'pts' using differences of forward and backward points
		/// (except for the first and last elements).
		/// </summary>
		/// <param name="pts"></param>
		/// <param name="output"></param>
		public static void CalcOrientation(PointF[] pts, double[] output, double kernelSd, bool wrap)
		{
			for(int i=0; i<pts.Length; i++)
				output[i] = KernelDerivative(pts, kernelSd, i, wrap);
		}
		public static double[,] ToDoubleArray(PointF[] pts)
		{
			double[,] output = new double[pts.Length,2];
			for(int i=0; i<pts.Length; i++)
			{
				output[i,0] = pts[i].X;
				output[i,1] = pts[i].Y;
			}
			return output;
		}
		public static int _end;
		/// <summary>
		/// Differences f into df using forward and backward values of f (even for first and last values of f, when 'wrap' is true).
		/// This assumes that f is an angular parameter which cannot go outside of [-pi, pi]
		/// </summary>
		/// <param name="f">input</param>
		/// <param name="df">output in [-pi, pi]</param>
		/// <param name="len">number of elements of 'f' to difference.</param>
		public static void AngularDiff(double[] f, double[] df, int len, bool wrap)
		{
			if(len>1)
			{
				_end = len-1;
				if(wrap)
				{
					//	First value
					df[0]=f[1]-f[_end];
					df[0]=Math.Atan2(Math.Sin(df[0]), Math.Cos(df[0]))/2.0;

					//	Last value
					df[_end]=f[0]-f[_end-1];
					df[_end]=Math.Atan2(Math.Sin(df[_end]), Math.Cos(df[_end]));
				}
				else
				{
					//	First value
					df[0]=f[1]-f[0];
					df[0]=Math.Atan2(Math.Sin(df[0]), Math.Cos(df[0]));

					//	Last value
					df[_end]=f[len-2]-f[_end];
					df[_end]=Math.Atan2(Math.Sin(df[_end]), Math.Cos(df[_end]));
				}
				//	Middle values.
				for(int i=1; i<len-2; i++)
				{
					df[i]=f[i+1]-f[i-1];
					df[i]=Math.Atan2(Math.Sin(df[i]), Math.Cos(df[i]))/2.0;
				}
			}
		}
		/// <summary>
		/// Derivative is returned as an angle.
		/// </summary>
		/// <param name="spine">PointF[]  A "spine" of adjacent points that make an edge or the center of a ribbon.</param>
		/// <param name="spineKernelSd">The standard deviation of the Gaussian smoothing kernel that smooths the derivative with respect
		/// to spine indices.</param>
		/// <param name="atIndex">The location on the spine where the derivative is to be evaluated.</param>
		/// <returns>The derivative returned is an angle (in radians).</returns>
		public static double KernelDerivative(PointF[] spine, double spineKernelSd, int atIndex, bool wrap)
		{
			int i;

			if(spine.Length<2)
				return 0.0;
			
			//	define the angle of the line
			if(spineKernelSd<=0.0f)	//	without kernel
			{
				i=atIndex;
				if(i==0)
					return Math.Atan2(spine[1].Y-spine[0].Y, spine[1].X-spine[0].X);
				else if(i==spine.Length-1)
					return Math.Atan2(spine[spine.Length-1].Y-spine[spine.Length-2].Y, spine[spine.Length-1].X-spine[spine.Length-2].X);
				else
					return Math.Atan2(spine[i+1].Y-spine[i-1].Y, spine[i+1].X-spine[i-1].X);
			}
			else //	with kernel
			{
				//	Total of weighted cos and sin values.
				double tCos, tSin;
				//	The current density
				double dens, z;
				//	The angle
				double theta;
				//	For wrapping
				int iBack, iFor;

				//	Start at the center and go backwards.
				i=atIndex;
				if(i==0)
					theta = Math.Atan2(spine[1].Y-spine[0].Y, spine[1].X-spine[0].X);
				else if(i==spine.Length-1)
					theta = Math.Atan2(spine[spine.Length-1].Y-spine[spine.Length-2].Y, spine[spine.Length-1].X-spine[spine.Length-2].X);
				else
					theta = Math.Atan2(spine[i+1].Y-spine[i-1].Y, spine[i+1].X-spine[i-1].X);
				dens = Util.Gaussian.Density( 0.0 );
				tCos=Math.Cos(theta)*dens;
				tSin=Math.Sin(theta)*dens;
				if(wrap)
				{
					while(true)
					{
						z = (double)(--i-atIndex)/spineKernelSd;
						if(z>-2.0)
						{
							iFor = i+1;
							iBack = i-1;
							while(iFor<0)
								iFor+= spine.Length;
							while(iBack<0)
								iBack+= spine.Length;
							dens = Util.Gaussian.Density( z );
							theta = Math.Atan2(spine[iFor].Y-spine[iBack].Y, spine[iFor].X-spine[iBack].X);
							tCos+=Math.Cos(theta)*dens;
							tSin+=Math.Sin(theta)*dens;
						}
						else
							break;
					}
				}
				else
				{
					while(--i>=0)
					{
						z = (double)(i-atIndex)/spineKernelSd;
						if(z > -2.0)	 //	only include within 2.0 sd
						{
							dens = Util.Gaussian.Density( z );
							if(i==0)
								theta = Math.Atan2(spine[1].Y-spine[0].Y, spine[1].X-spine[0].X);
							else
								theta = Math.Atan2(spine[i+1].Y-spine[i-1].Y, spine[i+1].X-spine[i-1].X);
							tCos+=Math.Cos(theta)*dens;
							tSin+=Math.Sin(theta)*dens;
						}
						else
							break;
					}
				}

				//	Go back to center and advance forwards.
				i=atIndex;
				if(wrap)
				{
					while(true)
					{
						z = (double)(++i-atIndex)/spineKernelSd;
						if(z<2.0)
						{
							iFor = i+1;
							iBack = i-1;
							while(iFor>=spine.Length)
								iFor-= spine.Length;
							while(iBack>=spine.Length)
								iBack-= spine.Length;
							dens = Util.Gaussian.Density( z );
							theta = Math.Atan2(spine[iFor].Y-spine[iBack].Y, spine[iFor].X-spine[iBack].X);
							tCos+=Math.Cos(theta)*dens;
							tSin+=Math.Sin(theta)*dens;
						}
						else
							break;
					}
				}
				else
				{
					while( ++i<spine.Length )
					{
						z = (double)(i-atIndex)/spineKernelSd;
						if(z < 2.0)	//	only include within 2.0 sd
						{
							dens = Util.Gaussian.Density( z );
							if(i==spine.Length-1)
								theta = Math.Atan2(spine[spine.Length-1].Y-spine[spine.Length-2].Y, spine[spine.Length-1].X-spine[spine.Length-2].X);
							else
								theta = Math.Atan2(spine[i+1].Y-spine[i-1].Y, spine[i+1].X-spine[i-1].X);
							tCos+=Math.Cos(theta)*dens;
							tSin+=Math.Sin(theta)*dens;
						}
						else
							break;
					}
				}
				return Math.Atan2(tSin, tCos);
			}
		}			
		public static float _dx1, _dx2, _dy1, _dy2, _dxTrail, _dyTrail, _dist;
		public static double _c1, _c2;
		public static PointF _trailer;
		/// <summary>
		/// Moves the points right-wards by the given distance.  If the last
		/// point is equal to the first point, then a closed polygon is assumed
		/// (you get inward movement if the polygon's order is standardized).
		/// </summary>
		/// <param name="pts">A list of points to be moved (this is replaced by the output).</param>
		/// <param name="moveDist">The right-ward (or standardized inward) distance to move.</param>
		public static void MoveRightwards(PointF[] pts, float moveDist)
		{
			//	Record the first point to later see if this is a closed polygon.
			_trailer = pts[0];

			if(pts.Length<=1)
				return;
			_dx2 = pts[1].X-pts[0].X;
			_dy2 = pts[1].Y-pts[0].Y;
			_dist = (float)Math.Sqrt(_dx2*_dx2+_dy2*_dy2);
			_dx2/=_dist;
			_dy2/=_dist;

			_dxTrail = _dx2;
			_dyTrail = _dy2;
			
			pts[0] = new PointF(
				pts[0].X-moveDist*_dy2,	//	dx is actually the negative of the change in y.
				pts[0].Y+moveDist*_dx2 );	//	dy is actually the change in x.

			_iLen = pts.Length-1;
			for(_i=1; _i<_iLen; _i++)
			{
				_dx1 = _dx2;
				_dy1 = _dy2;
				_dx2 = pts[_i+1].X-pts[_i].X;
				_dy2 = pts[_i+1].Y-pts[_i].Y;
				_dist = (float)Math.Sqrt(_dx2*_dx2+_dy2*_dy2);
				_dx2/=_dist;
				_dy2/=_dist;
				pts[_i] = new PointF(
					pts[_i].X-moveDist*0.5f*(_dy1+_dy2),	//	dx is actually the change in y.
					pts[_i].Y+moveDist*0.5f*(_dx1+_dx2) );	//	dy is actually the negative of the change in x.
			}
			if(_trailer==pts[_iLen])	//	closed polygon, recompute pts[0]
			{
				pts[_iLen] = pts[0] = new PointF(
					_trailer.X-moveDist*0.5f*(_dyTrail+_dy2),	//	dx is actually the change in y.
					_trailer.Y+moveDist*0.5f*(_dxTrail+_dx2) );	//	dy is actually the negative of the change in x.
			}
			else	//	straight edge.
			{
				pts[_iLen] = new PointF(
					pts[_iLen].X-moveDist*_dy2,	//	dx is actually the change in y.
					pts[_iLen].Y+moveDist*_dx2 );	//	dy is actually the negative of the change in x.
			}
		}

		public class DistanceF
		{
			public static float dx, dy;
			public static float Distance(PointF p1, PointF p2)
			{
				dx = p1.X-p2.X;
				dy = p1.Y-p2.Y;
				return (float)Math.Sqrt(dx*dx+dy*dy);
			}
		}
		public delegate void ProgressUpdater(float progress);
		public static float[] LinSpace(float low, float high, float increment)
		{
			int length = (int)(Math.Abs((high-low)/increment)+0.5f)+1;
			float[] output = new float[length];
			for(int i=0; i<length; i++)
				output[i] = low + (float)i*increment;
			return output;
		}
		public static float[] LinSpace(float low, float high, int length)
		{
			float increment = (high-low)/(float)(length-1);
			float[] output = new float[length];
			for(int i=0; i<length; i++)
				output[i] = low + (float)i*increment;
			return output;
		}
		/// <summary>
		/// Factors using the simple sieve approach.  Based on Matlab's factor.m but improved.
		/// </summary>
		/// <param name="n">The number to be factored.</param>
		/// <returns>A vector of factors.</returns>
		public static int[] Factor(int n)
		{
			if(n<4)
				return new int[] {n};

			int[] p = Primes(n);

			if(n==p[p.Length-1])
				return new int[] {p[p.Length-1]};

			bool found=true;
			int i;

			ArrayList output = new ArrayList();
			while(n>1 && found)
			{
				found = false;
				//	Find any prime numbers that are divisible.
				for(i=0; i<p.Length; i++)
				{
					if( n%p[i] == 0 )
					{
						found=true;
						output.Add(p[i]);
						n/=p[i--];
					}
				}
			}
			return (int[])output.ToArray(typeof(int));
		}
		/// <summary>
		/// Generates a vector of prime numbers which are less than or equal to n.
		/// Code is based on Matlab's primes.m.
		/// </summary>
		/// <param name="n">The greatest number that can possibly be returned.</param>
		/// <returns>A vector of prime numbers.</returns>
		public static int[] Primes(int n)
		{
			if( n<2 )
				return new int[0];

			int i,k;

			//	The number of positive odd numbers less than or equal to n.
			int q = (n+1)/2;

			//	The list of odd numbers, also '2' is the first prime number (instead of '1')
			int[] possibilities = new int[q];
			possibilities[0] = 2;
			for(i=1; i<q; i++)
				possibilities[i] = 2*i+1;

			//	The largest possible factor of n.
			int sqrtn = (int)Math.Sqrt((double)n);

			//	Some algorithm taken from matlab's primes.m
			for(k=3; k<=sqrtn; k+=2)
				if( possibilities[(k+1)/2-1] > 0 )
					for( i=(k*k+1)/2-1; i<q; i+=k)
						possibilities[i] = 0;

			//	All the non-zeros represent the answer.  Count them.
			int ct=0;
			for(i=0; i<q; i++)
				if(possibilities[i]>0)
					ct++;

			//	Generate output.
			int[] output = new int[ct];
			i=-1;
			k=0;
			while(k<ct)
				if(possibilities[++i]>0)
					output[k++]=possibilities[i];
			return output;
		}
		/// <summary>
		/// Computes the mean of each column of "data"
		/// </summary>
		/// <param name="data"></param>
		/// <returns>Returns a vector where each element is the mean for a data column.</returns>
		public static double[] Means(double[,] data)
		{
			int i,j;
			int rows = data.GetLength(0);
			int cols = data.GetLength(1);
			double[] output = new double[cols];
			for(i=0; i<cols; i++)
			{
				for(j=0; j<rows; j++)
				{
					output[i]+=data[j,i];
				}
				output[i]/=(double)rows;
			}
			return output;
		}
		/// <summary>
		/// Computes standard deviation of each column of "data" given the mean of each data column.
		/// </summary>
		/// <param name="data">The data</param>
		/// <param name="mean">The mean of each data column (i.e. the output of Means)</param>
		/// <returns>Returns a vector where each element is the standard deviation for a data column.</returns>
		public static double[] Sds(double[,] data, double[] mean)
		{
			int i,j;
			int rows = data.GetLength(0);
			int cols = data.GetLength(1);
			double[] output = new double[cols];
			double dx;
			if(rows > 1)
			{
				for(i=0; i<cols; i++)
				{
					for(j=0; j<rows; j++)
					{
						dx = (data[j,i]-mean[i]);
						output[i]+=dx*dx;
					}
					output[i]= Math.Sqrt( output[i] / (double)(rows-1) );
				}
			}
			return output;
		}
		/// <summary>
		/// Standardizes each column of data by subtracting its mean and dividing by its standard deviation.
		/// </summary>
		/// <param name="data">The data to be standardized</param>
		/// <param name="mean">The mean of each column</param>
		/// <param name="sd">The standard deviation of each column</param>
		public static void Standardize(double[,] data, double[] mean, double[] sd)
		{
			int i, j;
			int rows = data.GetLength(0);
			int cols = data.GetLength(1);
			for(i=0; i<cols; i++)
				for(j=0; j<rows; j++)
					data[j,i] = (data[j,i]-mean[i])/sd[i];
		}
		public static ArrayList Unique(ArrayList list)
		{
			if(list==null || list.Count==0)
				return null;
			list.Sort();
			bool[] unique = new bool[list.Count];
			unique[0] = true;
			int i;
			int ct=1;
			for(i=1; i<unique.Length; i++)
			{
				if( !(list[i]==list[i-1]) )
				{
					ct++;
					unique[i] = true;
				}
			}
			object[] output = new object[ct];
			while(ct>0)
			{
				if( unique[--i] )
					output[--ct] = list[i];
			}
			return new ArrayList(output);
		}
		/// <summary>
		/// Windows the values in io with the first derivative of a 2-d gaussian with the specified angle and standard deviation.
		/// </summary>
		/// <param name="sd">The standard deviation of the Gaussian.</param>
		/// <param name="angle">The angle of the Gaussian.</param>
		/// <param name="io">The input and output.</param>
		public static void GaussWindow1(double sd, double angle, double[,] io)
		{
			int ni=io.GetLength(0), nj=io.GetLength(1);
			double ic = (double)ni/2.0;
			double jc = (double)nj/2.0;
			double sin = Math.Sin(angle);
			double cos = Math.Cos(angle);
			double iy, jx, y, x, arg;
			double scale = 1.0/sd/sd;
			int i,j;
			for(i=0; i<ni; i++)
				for(j=0; j<nj; j++)
				{
					iy = (double)i-ic;
					jx = (double)j-jc;
					if(jx==0.0 && iy==0.0)
						io[i,j]=0.0;
					else
					{
                        //	Rotataion
                        y = iy*cos - jx*sin;
                        x = iy*sin + jx*cos;
                        //	Gradient
                        arg = y/sd;
						io[i,j] *= scale*arg*Math.Exp(-(x*x+y*y)*scale/2.0);
					}
				}
		}
		/// <summary>
		/// Windows the values in io with the second derivative of a 2-d gaussian with the specified angle and standard deviation.
		/// </summary>
		/// <param name="sd">The standard deviation of the Gaussian.</param>
		/// <param name="angle">The angle of the Gaussian.</param>
		/// <param name="io">The input and output.</param>
		public static void GaussWindow2(double sd, double angle, double[,] io)
		{
			int ni=io.GetLength(0), nj=io.GetLength(1);
			double ic = (double)ni/2.0;
			double jc = (double)nj/2.0;
			double sin = Math.Sin(angle);
			double cos = Math.Cos(angle);
			double iy, jx, y, x, arg;
			double scale = 1.0/sd/sd;
			int i,j;
			for(i=0; i<ni; i++)
				for(j=0; j<nj; j++)
				{
					iy = (double)i-ic;
					jx = (double)j-jc;
					if(jx==0.0 && iy==0.0)
						io[i,j]=0.0;
					else
					{
                        //	Rotataion
                        y = iy*cos - jx*sin;
                        x = iy*sin + jx*cos;
                        //	Gradient
                        arg = y/sd;
                        io[i,j] *= scale*(arg*arg-scale)*Math.Exp(-(x*x+y*y)*scale/2.0);
					}
				}
		}
        /// <summary>
        /// Sets the elements of an array to a specified value.  This uses 10:1 "loop unrolling" for faster performance.
        /// </summary>
        /// <param name="array"></param>
        /// <param name="val"></param>
        public static void SetElements(float[,] array, float val)
        {
            int ni = array.GetLength(0);
            //int nj = array.GetLength(1)/10;
            int nj = array.GetLength(1);
            int jstride=0;
            int i, j;
            for (i = 0; i < ni; i ++)
            {
                for (j = 0; j < nj; j++)
                    array[i, j] = val;
                /*
                for (j = 0; j < nj; j ++)
                {
                    jstride = j * 10;
                    array[i, jstride] = val;
                    array[i, jstride + 1] = val;
                    array[i, jstride + 2] = val;
                    array[i, jstride + 3] = val;
                    array[i, jstride + 4] = val;
                    array[i, jstride + 5] = val;
                    array[i, jstride + 6] = val;
                    array[i, jstride + 7] = val;
                    array[i, jstride + 8] = val;
                    array[i, jstride + 9] = val;
                }
                jstride = nj*10;
                nj = array.GetLength(1)%10;
                for (j = 0; j < nj; j++)
                    array[i, jstride + j] = val;
                */
            }
        }
    }
	/// <summary>
	/// A minimization function written by me for a certain class of problems that Numerical Recipes doesn't solve quickly.
	/// This is a relatively stupid solver (no Jacobian) with an adaptive step-size.
	/// </summary>
	public class IngMin
	{
		public string FullVerbose
		{
			get
			{
				int i;
				string output = ChisqSol.ToString() + " +/- " + ChisqErrTol.ToString() + '\n' + Ngoods.ToString() + " good steps taken\n";
				for(i=0; i<A.Length; i++)
					output = output + A[i].ToString() + " --> " + Asol[i].ToString() + '\n';
				return output;
			}
		}
		public int dim;
		public double chisq, chisqTry, firstScalar, scalar, fullScale, low, high, tol, fac;
		public double[] a, aTry, aStep, aStepLast, temp;
		public delegate double DelChisq(double[] a);
		public DelChisq Func;
		public double ChisqErrTol;
		public double[] A, Asol;
		public double ChisqSol;
		public int Ngoods;
		public IngMin(double ChisqErrTol, double[] A, DelChisq Func)
		{
			this.ChisqErrTol = ChisqErrTol;
			dim = A.Length;
			this.A = A;
			a = new double[dim];
			A.CopyTo(a,0);
			aTry = new double[dim];
			aStep = new double[dim];
			aStepLast = new double[dim];
			for(int i=0; i<dim; i++)
				aStep[i] = aStepLast[i] = 1.0;
			tol = ChisqErrTol/(double)dim;
			this.Func = Func;
		}
		public int Solve()
		{
			ChisqSol = chisq = Func(a);
			int nGood = 0;
			int slowGood=1, turboGood=1;
			while(slowGood+turboGood > 0)
			{
				slowGood = SlowStep(0.1);
				turboGood = TurboStep();
				nGood+=slowGood+turboGood;
			}
			Asol = new double[dim];
			a.CopyTo(Asol,0);
			Ngoods = nGood;
			return nGood;
		}
		/// <summary>
		/// Turbo minimization mode.  Keep taking aStep * scalar until there is no more minimization.
		/// The member 'a' is set to the turbo solution, aStepLast is set to the last turbo step size (aStep * scalar * fullStepSize)
		/// where fullStepSize is the scale of aStep that was taken.
		/// </summary>
		/// <returns>The number of successful function evals (which does not necessarily scale exactly with the change in a).</returns>
		public int TurboStep()
		{
			int i, nGood=0;
			chisq = ChisqSol;
			firstScalar = scalar = 2.0;
			for(i=0; i<dim; i++)
			{
				scalar = firstScalar;
				while(scalar>=firstScalar)
				{
					aTry[i] = a[i] + aStep[i]*scalar;
					chisqTry = Func(aTry);
					if(chisqTry<chisq)
					{
						if(chisqTry+tol<chisq)
						{
							nGood++;
							aStepLast[i] = aStep[i]*scalar;
							scalar*=scalar;
						}
						else
						{
							scalar = Math.Sqrt(scalar);
						}

						a[i] = aTry[i];
						chisq = chisqTry;
					}
					else
					{
						scalar = Math.Sqrt(scalar);
					}
				}
				aTry[i] = a[i];
			}
			ChisqSol = chisq;
			return nGood;
		}
		public int SlowStep(double minRangeFac)
		{
			int nGood = 0;
			a.CopyTo(aTry,0);
			int i;
			for(i=0; i<dim; i++)
			{
				fac = 1.0;
				while(fac>minRangeFac)
				{
					aStepLast[i] = -aStepLast[i];
					aTry[i] = a[i] + aStepLast[i]*fac;
					chisqTry = Func(aTry);
					if(chisqTry<chisq)
					{
						if(chisqTry+tol<chisq)
							nGood++;
						a[i] = aTry[i];
						aStep[i] = aStepLast[i];
						chisq = chisqTry;
						ChisqSol = chisq;
						fac = 0.0;
					}
					else
					{
						fac/=2.0;
						aStepLast[i] = -aStepLast[i];
						aTry[i] = a[i] + aStepLast[i];
						chisqTry = Func(aTry);
						if(chisqTry<chisq)
						{
							if(chisqTry+tol<chisq)
								nGood++;
							a[i] = aTry[i];
							aStep[i] = aStepLast[i];
							chisq = chisqTry;
							ChisqSol = chisq;
							fac = 0.0;
						}
						else
						{
							fac/=4.0;
							aTry[i] = a[i];
						}
					}
				}
				if(fac>0.0)
				{
					aStep[i] *= fac;
					aStepLast[i] = aStep[i];
				}
			}
			return nGood;
		}
	}
}


